06_analyze_DEG_and_signatures
Amy Huang
2024-03-07
Read in data and results
gene_count_D9_df_filtered <- readRDS(here("bulkRNAseq", "results", "gene_count_D9_df_filtered.rds"))
gene_count_D30_df_filtered <- readRDS(here("bulkRNAseq", "results", "gene_count_D30_df_filtered.rds"))
sample_D9_metadata_df <- read_tsv(here("bulkRNAseq", "results", "sample_metadata_D9_df.txt"))## Rows: 23 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): sample, day, mice, genotype, exh_subset, sample_num
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.sample_D30_metadata_df <- read_tsv(here("bulkRNAseq", "results", "sample_metadata_D30_df.txt"))## Rows: 12 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): sample, day, mice, genotype, exh_subset, sample_num
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.norm_counts_D9 <- readRDS(here("bulkRNAseq", "results", "deseq_norm_counts_D9_exh.prog.vs.term.rds")) %>%
as_tibble(rownames = "gene")
norm_counts_D30 <- readRDS(here("bulkRNAseq", "results", "deseq_norm_counts_D30_exh.prog.vs.term.rds")) %>%
as_tibble(rownames = "gene")dge_D30_compare_exh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_exh.prog.vs.term.txt")) # WT## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D9_compare_exh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_exh.prog.vs.term.txt")) # WT## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D9_progexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_progexh_genotype.Enhdel.v.WT.txt"))## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D9_transexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_transexh_genotype.Enhdel.v.WT.txt"))## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D9_termexh <- read_tsv(here("bulkRNAseq", "results", "dge_D9_termexh_genotype.Enhdel.v.WT.txt"))## Rows: 23255 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D30_progexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_progexh_genotype.Enhdel.v.WT.txt"))## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D30_transexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_transexh_genotype.Enhdel.v.WT.txt"))## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.dge_D30_termexh <- read_tsv(here("bulkRNAseq", "results", "dge_D30_termexh_genotype.Enhdel.v.WT.txt"))## Rows: 18983 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gene
## dbl (6): baseMean, log2FoldChange, lfcSE, stat, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.sc_dge <- read_tsv(here("scRNAseq", "data", "processed_data_tables", "01_dge.genotype_one-v-one_per-cluster.tsv"))## Rows: 216930 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): gene, contrast, Clusters
## dbl (5): p_val, avg_logFC, pct.1, pct.2, p_val_adj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.gs_df <- read_tsv(here("scRNAseq", "data", "gene_signatures",
"02_gene_signatures.msigdbr_H-C2-C7.tsv"))## Rows: 1509466 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (13): gs_cat, gs_subcat, gs_name, gene_symbol, ensembl_gene, human_gene_...
## dbl (5): entrez_gene, human_entrez_gene, gs_pmid, taxon_id, num_ortholog_so...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.gs_subset_bp_df <- gs_df %>%
filter(str_detect(gs_name, "HALLMARK|REACTOME"))
length(unique(gs_subset_bp_df$gs_name))## [1] 1665# 1665
gs_subset_bp_list <- gs_subset_bp_df %>%
named_group_split(gs_name) %>%
map(~..1$gene_symbol)
gs_subset_tcell_df <- gs_df %>%
filter(str_detect(gs_name, "_T_CELL|TCELL|CD8|CD4|EXH")) %>%
filter(str_detect(gs_name, "HALLMARK|REACTOME", negate = T)) %>%
filter(str_detect(gs_name, "NEUTROPHIL|MONO|ERYTHROBLAST|DC|NKT|LIN|BCELL|LSK|CD40|CD44|MAST", negate = T))
length(unique(gs_subset_tcell_df$gs_name))## [1] 1446# 1446
gs_subset_tcell_list <- gs_subset_tcell_df %>%
named_group_split(gs_name) %>%
map(~..1$gene_symbol)gsea_D9_tcell <- read_tsv(
here("bulkRNAseq", "results", "gsea_tcell_wide_D9_genotype.Enhdel.v.WT.txt")
)## Rows: 1445 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gs_name
## dbl (12): p.val_progexh, p.val_termexh, p.val_transexh, q.val_progexh, q.val...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.gsea_D9_bp <- read_tsv(
here("bulkRNAseq", "results", "gsea_bp_wide_D9_genotype.Enhdel.v.WT.txt")
)## Rows: 1660 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): gs_name
## dbl (12): p.val_progexh, p.val_termexh, p.val_transexh, q.val_progexh, q.val...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.DGE, term vs prog exh
my_plot_volcano <- function(dge_df, label_genes) {
dge_df %>%
drop_na(padj) %>%
ggplot() +
aes(log2FoldChange, -log10(padj)) +
geom_point(size = 0.5, alpha = 0.5) +
geom_hline(yintercept = -log10(0.05), linetype = "dotted", color = "red") +
geom_point(aes(color = padj < 0.05),
data = ~..1 %>% filter(gene %in% c(label_genes)),
size = 2) +
geom_text_repel(aes(label = gene, color = padj < 0.05),
data = ~..1 %>% filter(gene %in% c(label_genes)),
max.overlaps = Inf) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
theme(legend.position = "none")
}label_genes <- c("Havcr2", "Slamf6")
dge_D30_compare_exh %>%
my_plot_volcano(label_genes = label_genes) +
labs(title= "D30")
dge_D9_compare_exh %>%
my_plot_volcano(label_genes = label_genes) +
labs(title= "D9")
# gene list from miller et al 2019
inh_receptors <- c("Pdcd1", "Havcr2", "Lag3", "Cd244", "Entpd1", "Cd38", "Cd101", "Tigit", "Ctla4")
surface_receptors <- c("Tnfsf9", "Tnfrsf4", "Il2rb", "Klrg1", "Cd28", "Icos", "Tnfsf14", "Sell", "Il7r")
eff_molecules <- c("Ifng", "Il10", "Gzma", "Gzmb", "Prf1", "Tnfsf10", "Fasl", "Tnf", "Il2")
tfs <- c("Nfkb2", "Id2", "Hif1a", "Bach2", "Nfkb1", "Id3", "Tcf7", "Lef1", "Satb1", "Batf",
"Nfatc1", "Prdm1", "Runx1", "Runx3", "Tbx21", "Tox", "Bcl6", "Eomes")
chemokines <- c("Cx3cr1", "Ccl5", "Ccl4", "Ccl3", "Csf1", "Cxcr5", "Ccr7", "Xcl1", "Cxcl10")
select_genes <- c(inh_receptors, surface_receptors, eff_molecules, tfs, chemokines)plot_gene_heatmap <- function(norm_counts, gene_list,
col_index_ordering = c(3, 6, 9, 11, 2, 5, 8, 10, 1, 4, 7),
hheight = unit(1.5, "in"),
hwidth = unit(2, "in")) {
h_mtx <-
norm_counts %>%
filter(gene %in% gene_list) %>%
column_to_rownames("gene") %>%
as.matrix() %>%
t() %>%
scale() %>%
t() %>%
.[rownames(.)[order(match(rownames(.), gene_list))],
col_index_ordering]
Heatmap(
h_mtx,
width = hwidth,
height = hheight,
col = colorRamp2(colors = c("purple4","yellow"),
breaks = c(-3, 3)),
name = "row scaled\nnorm counts",
cluster_rows = FALSE,
cluster_columns = FALSE
)
}plot_gene_heatmap(norm_counts_D9, inh_receptors)
plot_gene_heatmap(norm_counts_D9, surface_receptors)
plot_gene_heatmap(norm_counts_D9, eff_molecules)
plot_gene_heatmap(norm_counts_D9, tfs[1:9])
plot_gene_heatmap(norm_counts_D9, tfs[10:18])
plot_gene_heatmap(norm_counts_D9, chemokines)
plot_gene_heatmap(norm_counts_D30, inh_receptors, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
plot_gene_heatmap(norm_counts_D30, surface_receptors, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
plot_gene_heatmap(norm_counts_D30, eff_molecules, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
plot_gene_heatmap(norm_counts_D30, tfs[1:9], c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
plot_gene_heatmap(norm_counts_D30, tfs[10:18], c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
plot_gene_heatmap(norm_counts_D30, chemokines, c(3, 6, 2, 5, 1, 4), hwidth = unit(1, "in"))
DGE bw genotypes
select_genes <- c("Btg2", "Ctse", "Cxcr4", "Il17ra", "Tnfrsf9")
dge_D9_progexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, prog exh")
dge_D9_transexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, trans exh")
dge_D9_termexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, term exh")
From single cell too:
select_genes <- c("Btg2", "Ctse", "Cxcr4", "Il17ra", "Tnfrsf9")
sc_dge %>%
drop_na(p_val_adj) %>%
filter(contrast == "D_W") %>%
filter(Clusters %in% c("terminal", "progenitor", "transitory")) %>%
ggplot() +
aes(avg_logFC,-log10(p_val_adj)) +
geom_point(size = 0.5, alpha = 0.5) +
geom_hline(
yintercept = -log10(0.05),
linetype = "dotted",
color = "red"
) +
geom_point(aes(color = p_val_adj < 0.05),
data = ~ ..1 %>% filter(gene %in% c(select_genes)),
size = 2) +
geom_text_repel(
aes(label = gene, color = p_val_adj < 0.05),
data = ~ ..1 %>% filter(gene %in% c(select_genes)),
max.overlaps = Inf
) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
theme(legend.position = "none") +
facet_grid(~Clusters) +
theme(strip.background = element_blank()) +
labs(title = "DGE from scRNAseq")
select_genes <- c("Sell", "Tcf7")
dge_D9_progexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, prog exh")
dge_D9_transexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, trans exh")
dge_D9_termexh %>% my_plot_volcano(label_genes = select_genes) +
labs(title = "D9, term exh")
From single cell too:
select_genes <- c("Sell", "Tcf7")
sc_dge %>%
drop_na(p_val_adj) %>%
filter(contrast == "D_W") %>%
filter(Clusters %in% c("terminal", "progenitor", "transitory")) %>%
ggplot() +
aes(avg_logFC,-log10(p_val_adj)) +
geom_point(size = 0.5, alpha = 0.5) +
geom_hline(
yintercept = -log10(0.05),
linetype = "dotted",
color = "red"
) +
geom_point(aes(color = p_val_adj < 0.05),
data = ~ ..1 %>% filter(gene %in% c(select_genes)),
size = 2) +
geom_text_repel(
aes(label = gene, color = p_val_adj < 0.05),
data = ~ ..1 %>% filter(gene %in% c(select_genes)),
max.overlaps = Inf
) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
theme(legend.position = "none") +
facet_grid(~Clusters) +
theme(strip.background = element_blank()) +
labs(title = "DGE from scRNAseq")
GSEA
# gsea_res_bp %>%
# dplyr::select(gs_name, str_subset(colnames(.), "prog")) %>%
# filter(q.val_progexh < 0.05) %>%
# arrange(desc(sscore_progexh)) %>%
# View()
# gsea_res_bp %>%
# dplyr::select(gs_name, str_subset(colnames(.), "trans")) %>%
# filter(q.val_transexh < 0.05) %>%
# arrange(sscore_transexh) %>%
# View()dge_list <- list(
progenitor = dge_D9_progexh %>% mutate(cluster_id = "progexh"),
transitory = dge_D9_transexh %>% mutate(cluster_id = "transexh"),
terminal = dge_D9_termexh %>% mutate(cluster_id = "termexh")
)
dge_df <- dge_list %>%
purrr::reduce(bind_rows) %>%
mutate(signed_neglog10_p = ifelse(log2FoldChange > 0, -log10(pvalue), log10(pvalue))) %>%
drop_na(log2FoldChange, pvalue)
make_rank_list <- function(.data,
gene = gene,
rank_by = log2FoldChange) {
rank_by <- enquo(rank_by)
gene <- enquo(gene)
data_fc <- .data %>%
arrange(desc(!!rank_by))
data_fc %>%
.[[rlang::as_name(rank_by)]] %>%
set_names(data_fc[[rlang::as_name(gene)]])
}
rank_lists <- dge_df %>%
named_group_split(cluster_id) %>%
map(~make_rank_list(..1, rank_by = signed_neglog10_p, gene = gene))my_plot_gsea_curve <- function(gsea_res, rank_list, gs_subset_list, gs_name_str, exh_cluster) {
my_gs <- gsea_res %>%
filter(gs_name == gs_name_str)
Rsc::get_leading_edge_df(rank_vec = rank_list[[exh_cluster]],
gs_vec = gs_subset_list[[gs_name_str]]) %>%
Rsc::plot_gsea_curve(title = paste0(gs_name_str, ", ", exh_cluster)) +
theme(legend.position = "none") +
annotate(
"text",
x = Inf,
y = Inf,
hjust = 1,
vjust = 1,
label = paste0(
"qval=",
format(my_gs[[paste0("q.val_", exh_cluster)]], digits = 2),
"\n",
"sscore=",
format(my_gs[[paste0("sscore_", exh_cluster)]], digits = 4)
)
)
}my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "progexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "transexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "REACTOME_CELL_CYCLE", "termexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "progexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "transexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_bp, rank_lists, gs_subset_bp_list, "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "termexh") +
ylim(-0.3, 0.3)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "progexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "progexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "transexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "transexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_DN", "termexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EXHAUSTED_VS_MEMORY_CD8_TCELL_UP", "termexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "progexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "progexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "transexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "transexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_DN", "termexh") +
ylim(-0.23, 0.23)
my_plot_gsea_curve(gsea_D9_tcell, rank_lists, gs_subset_tcell_list, "GSE9650_EFFECTOR_VS_EXHAUSTED_CD8_TCELL_UP", "termexh") +
ylim(-0.23, 0.23)
gene_sets <- str_subset(names(gs_subset_bp_list), "RUNX")
gsea_curve_list <- list()
for (gs in gene_sets) {
gsea_curve_list[[paste0("progexh / ", gs)]] <-
my_plot_gsea_curve(gsea_D9_bp, rank_lists,
gs_subset_bp_list, gs, "progexh")
gsea_curve_list[[paste0("termexh / ", gs)]] <-
my_plot_gsea_curve(gsea_D9_bp, rank_lists,
gs_subset_bp_list, gs, "termexh")
}gsea_curve_list[c("progexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY",
"termexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY",
"progexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY",
"termexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY",
"progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2",
"termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2",
"progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1",
"termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1")]## $`progexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY`
##
## $`termexh / REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY`
##
## $`progexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY`
##
## $`termexh / REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY`
##
## $`progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2`
##
## $`termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX2`
##
## $`progexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1`
##
## $`termexh / REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1`
Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
gs_vec = gs_subset_bp_list[["REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1"]]) %>%
# filter(is_leading_edge) %>%
filter(in_gs) %>%
write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_TRANSCRIPTIONAL_REGULATION_BY_RUNX1.txt"))Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
gs_vec = gs_subset_bp_list[["REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY"]]) %>%
# filter(is_leading_edge) %>%
filter(in_gs) %>%
write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_REGULATION_OF_RUNX2_EXPRESSION_AND_ACTIVITY.txt"))Rsc::get_leading_edge_df(rank_vec = rank_lists[["progexh"]],
gs_vec = gs_subset_bp_list[["REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY"]]) %>%
# filter(is_leading_edge) %>%
filter(in_gs) %>%
write_tsv(here("bulkRNAseq", "results", "leadingedge_progexh_D9_REACTOME_REGULATION_OF_RUNX3_EXPRESSION_AND_ACTIVITY.txt"))Rsc::get_leading_edge_df(rank_vec = rank_lists[["termexh"]],
gs_vec = gs_subset_tcell_list[["GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"]]) %>%
filter(is_leading_edge) %>%
filter(in_gs)## # A tibble: 71 × 7
## name value in_gs rank running_es direction is_leading_edge
## <chr> <dbl> <lgl> <int> <dbl> <chr> <lgl>
## 1 Dusp2 -0.501 TRUE 18517 -0.227 dn TRUE
## 2 H2ax -0.512 TRUE 18586 -0.226 dn TRUE
## 3 Anln -0.518 TRUE 18632 -0.223 dn TRUE
## 4 Nusap1 -0.542 TRUE 18757 -0.223 dn TRUE
## 5 Tmem97 -0.557 TRUE 18824 -0.222 dn TRUE
## 6 Ak3 -0.558 TRUE 18830 -0.217 dn TRUE
## 7 Dtl -0.560 TRUE 18839 -0.212 dn TRUE
## 8 Tuba1a -0.562 TRUE 18846 -0.207 dn TRUE
## 9 Plk4 -0.566 TRUE 18870 -0.203 dn TRUE
## 10 Nudt4 -0.567 TRUE 18878 -0.199 dn TRUE
## # … with 61 more rowsgs <- "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"
my_plot_gsea_curve(gsea_D9_tcell, rank_lists,
gs_subset_tcell_list, gs, "progexh")
my_plot_gsea_curve(gsea_D9_tcell, rank_lists,
gs_subset_tcell_list, gs, "termexh")